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Abstract. Brain rhythms contribute to every aspect of brain function. Here, we study critical and 
resonance phenomena that precede the emergence of brain rhythms. Using an analytical approach 
and simulations of a cortical circuit model of neural networks with stochastic neurons in the pres- 
ence of noise, we show that spontaneous appearance of network oscillations occurs as a dynamical 
(non-equilibrium) phase transition at a critical point determined by the noise level, network struc- 
ture, the balance between excitatory and inhibitory neurons, and other parameters. We find that the 
relaxation time of neural activity to a steady state, response to periodic stimuli at the frequency of 
the oscillations, amplitude of damped oscillations, and stochastic fluctuations of neural activity are 
dramatically increased when approaching the critical point of the transition. 
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INTRODUCTION 

Brain rhythms contribute in every aspect of brain function from sensory and cognitive 
processing, and memory to motor control [1]. Origin and physiological functions of 
brain rhythms are a topic problem in neuroscience. Brain rhythms are also related to 
many unusual phenomena observed in the brain. Interactions between billions of neurons 
give rise to phase transitions, self-organization, and critical phenomena [2, 3]. Phase 
transitions were observed, for example, in human bimanual coordination [4, 5, 6, 7, 8] 
and in living neural networks stimulated by electric fields [9]. There are evidences that 
epileptic seizures, alpha and gamma oscillations, and the ultraslow oscillations of BOLD 
fMRI patterns emerge as a result of non-equilibrium phase transitions. Neural avalanches 
are one more example of critical collective phenomena observed in the brain [10, 3]. 

Various resonance phenomena were also observed in the brain. Experimental inves- 
tigations of CA1 neuronal networks from mammalian brain demonstrated that stochas- 
tic resonance can enhance effects of intrinsic 4-10 Hz hippocampal theta and 40 Hz 
gamma oscillations [11]. Recently, using a functional imaging technique, Sasaki et al. 
[12] revealed that the majority of rat CA1 neurons act collectively like a band-pass filter. 
Damped oscillations and the Berger effect are also related to brain rhythms. The Berger 
effect manifests itself in activation of alpha waves on the electroencephalogram when 
the eyes are closed and diminution of alpha waves when they are opened [13]. 

In the present paper, we study collective dynamics of neural networks composed by 
excitatory and inhibitory neurons in the presence of noise. Based on exact analytical cal- 
culations and numerical simulations, we show that spontaneous emergence of network 



oscillations occurs as a dynamical (non-equilibrium) phase transition at a critical level 
of noise. The transition manifests itself in slowing down of the relaxation of a perturbed 
neural activity to a steady state, a strong enhancement of stochastic fluctuations of ac- 
tivities of neural populations and an increase of the linear response function to afferent 
periodic stimuli at the frequency of neural oscillations. We show that near to the critical 
boundary, neural networks act as damped harmonic oscillators or band-pass filters that 
pass frequencies within a certain range and attenuate frequencies outside that range. 



CORTICAL CIRCUIT MODEL 

We use a cortical circuit model [14] composed of N e pyramidal cells (excitatory neurons) 
and A 7 , interneurons (inhibitory neurons) that form a sparsely connected network. The 
probability that there is a synaptic connection between two neurons is c/N where 
N = N e +Ni is the total number of neurons and c is the mean degree. This network has the 
structure of a directed classical random graph (or Erdos-Renyi graph) with the Poisson 
degree distribution P n (c) = c n e~ c ' jn\ where n is the number of presynaptic neurons. 
Neurons receive sporadic inputs from a remote part of the cortex and synaptic noise. 
Neurons fire with a constant firing frequency v that is the same for both excitatory and 
inhibitory neurons. The total input V m to a neuron with index m, m = 1,2, . . . ,N, is the 
sum of random spikes from noise, excitatory and inhibitory neurons, 

N 

V m(t) = Y, k n( t ) a nmJnm + ^(t), (1) 
n=l 

where k n (t) is the number of spikes that arrive from presynaptic neuron n during the 
time interval [t — x,t], % is the integration time. Below we will consider the case TV < 1 
when the number of spikes k n (t) is 1 or 0. If we assume that the emissions times of 
spikes of different neurons are uncorrelated, then the parameter TV has a meaning of 
the probability that a postsynaptic neuron receives a spike from an active presynaptic 
neuron during time T. Furthermore, a nm is the adjacency matrix, i.e., a nm = 1 if there is 
a direct edge from neuron n to neuron m, otherwise a nm = 0. J nm is the efficacy of the 
synapse connecting neuron n with neuron m. J nm is positive if presynaptic neuron n is 
excitatory and it is negative if the neuron is inhibitory. £, (t) is the number of random 
spikes from noise that neuron m receives during the time interval [t — x,t\. We use the 
Gaussian distribution for ^ (t), 



G(|) =Aexp 



2o 2 
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where A is the normalization constant, o 2 is the variance, (n) is the mean number of 
random spikes determined by the mean rate CO rs , in) = (O rs x. Note that noise in our 
model is actually shot noise. According to Schottky's theorem, the intensity of this noise 
is proportional to (n). 



We consider stochastic neurons. Their response on input is a stochastic process that 
occurs with a certain rate. Two rules determine dynamics of stochastic neurons [14]: 

1. If the total input V m {t) at an inactive excitatory or inhibitory neuron m at time t is 
at least a certain threshold Q. (i.e., V m (t) > Q.), then this neuron is activated at a rate 
\x e or fa, respectively, and fires with a cyclic frequency v. 

2. Active excitatory (inhibitory) neuron m is inactivated at a rate fl e (fid if V m (t) < £1. 

We assume that l/p e and are of the order of the first spike latencies of excitatory 
and inhibitory neurons, respectively. We introduce the ratio 

a = Hi/ n e (3) 

that plays an important role in our model, as it will be shown below. The advantage of 
this model with stochastic neurons is that it can be solved analytically. 

In numerical simulations, we studied sparsely connected networks of size N = 10 3 — 
10 5 and applied the following algorithm. We divided time t into intervals of width 
At = T. At each time step, for each neuron, we calculated the input Eq. (1), taking into 
account that each active presynaptic neuron contributes with a spike with probability 
TV. The number of random spikes from noise in this input is generated according to the 
Gaussian distribution, Eq. (2). Then, with the probability T/i a , a = e,i, we updated the 
states of all neurons using the stochastic rules formulated above. We used the following 
parameters: the fraction of excitatory neurons is g e = N e /N = 75%, the fraction of 
inhibitory neurons is g, = Nj/N = 25%, the mean number of connections c = 1000 (750 
excitatory and 250 inhibitory connections), the threshold Q. = 30, and the variance of 



noise o 2 = 10. Following [15], we chose 7 !e = Ju = 7,, J ee = J e i = J e , and 7, = — 37 ( 



These parameters agree with anatomical estimates for cortex. In cortex, the fraction gi 
of inhibitory neurons is between 0.15 and 0.3, the mean number of synaptic connections 
c is about 7000. The threshold Q. is between 15 and 30 in neural networks in vivo [9] 
and about 30 — 400 in the brain. The level of noise in) was varied in the interval 0—150 
spikes per integration time x. We also assumed that, for simplicity, TV = 1 and T/i e = 0.1. 

Dynamical behavior of the model is described by the fractions p e {t) and p,-(f ) of active 
excitatory and inhibitory neurons, respectively, at time t. We will call them 'activities' 
of the neural populations. Using the rules of the stochastic dynamics formulated above 
and assuming that activities are changed slightly during the integration time T, in the 
infinite size limit N — > °°, we find a rate equation [14], 

d ^-=fa(t)(l- Pa (t))-p a (t)+^ a (p e (t),p i (t)). (4) 
jladt 

for a = e,i. The function m a (p e ,pi) is the probability that at time t the input to a 
randomly chosen excitatory or inhibitory neuron is at least the threshold £1. For the 
model under consideration ^-(p^p;) = m e {p e ,Pi) = ^(pe^Pi), where 

CO CO oo 

%,Pi) = II £ ®(Jek+JilH-0)G^)Pk(8ePeC)Pl(8iPiC). (5) 

fc=0Z=0|=-oo 

Here ®{x) is the Heaviside step function, the parameter c is defined as c = cvt, and 
Pk(gepec) and Pi(gipic) are the probabilities that a randomly chosen neuron receives k 
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FIGURE 1. Schematic phase diagram of the cortical model and critical and resonance phenomena near 
the critical boundary of the non-equilibrium phase transition to sustained network oscillations. 



spikes from active presynaptic excitatory and I spikes from inhibitory neurons, respec- 
tively, during the time window x at given activities p e and p,. The functions f e (t) and 
fi(t) represent a rate of spontaneous activation of excitatory and inhibitory neurons, re- 
spectively, by stimulus, for example, an electric field. The rate equation (4) is similar to 
the Wilson-Cowan equations [16, 17], see also [14]. Equation (4) is asymptotically exact 
in the limit N — > °°. 

Steady states of the neural populations can be found from Eq. (4), supposing dp a /dt = 
in the limit t — > °°. If p e (t) and p,-(f) at time t are close to steady state activities p e (°°) 
and p,(°°), then Eq. (4) enables us to describe relaxation of p a (t) to the steady state. We 
introduce 

5p a (t) = p a {t) - p a H = Re(A a e~' Yt ) (6) 

where A a is a complex amplitude. Using the standard perturbation theory, we solve 
Eq. (4) in the first order in 8p a (t). We find 



y± = +B 2 ) ± l - (B { - B 2 r AaD cl D lt 



1/2 

(7) 



where we introduced parameters B\ = 1 —D ee , B 2 = Ct(l —Da), D a t, = d x ¥ a (p e: pi)/dpt, 
for a,b = e,i. respectively. Derivatives D a b are determined by the activities p e (°°) and 
Pi(°°) from the non-linear equation Eq. (4) when dp a /dt = [14]. The real and imagi- 
nary parts of the complex rate y (y r = Re(j-) and y ; = Im{yJ)) determine the relaxation 
rate and the angular frequency of damped oscillations, respectively. Notice that the pe- 
riod of the oscillations equals iTz/ji. 

Analyzing behavior of y r and y, in dependence on a and (n), we obtain the phase 
diagram in Fig. 1. One can see that there are three regions. There is a region I (small 
noise level and/or large a) where the relaxation of the neural activity to a steady state 
is exponential (y r > and y- = 0). In region II, the neural activity relaxes in a form of 
damped oscillations (y r > and y ■ ^ 0). In region III, network oscillations are sustained. 



A similar phase diagram was found in [14] for a simpler model. If a is above a critical 
value a t , that corresponds to the a-coordinate of the top point of the region III in Fig. 1, 
then with increasing the noise level (n), the activities p e and p, in the steady state 
undergo a first-order phase transition at a critical noise level n c . A similar discontinuous 
transition was observed in living neural networks in vitro when living neural networks 
were stimulated by an electric field [9]. Neuronal avalanches are precursors of this phase 
transition. Activation (or inactivation) of one neuron can trigger avalanche process of 
activation (or inactivation) of a cluster of neurons. In cortex, neuronal avalanches have 
been observed experimentally [10], see the review [3]. 

If the parameter a < a t , sustained networks oscillations appear in a certain 'optimal' 
range of the noise level (n) between two critical points. Weak noise can not stimulate 
network oscillations. Too strong noise over-activates neural networks and only damped 
oscillations can occur. The critical boundary of region with the sustained oscillations is 
determined by the condition that the relaxation rate y r is zero, 

Yr = Re(r-)=0, (8) 

where the complex frequency y_ is given by Eq. (7). For the parameters given above and 
T = 10 ms, frequencies of the oscillations lie in the range of brain waves (1- 100 Hz). 



LINEAR RESPONSE FUNCTION AND BAND-PASS FILTER 

BEHAVIOR 

Now we study critical phenomena that precede the non-equilibrium phase transition 
from asynchronious dynamics to sustained oscillations. For this purpose we calculate 
the linear response of the neural network to a time-dependent stimulus f e (t) and fi(t) in 
Eq. (4) for region I and II. Here we are not studying a response in region III that needs 
a special consideration. A response of the neural population a = e, i to a weak stimulus 
f a (t) is determined by the linear response function %ab{t — t'), 

Ap a (t)=p a (t)-p a (oo)=^ f Xab(t-t')f b (t')dt'. (9) 

b=ef- oa 

Solving Eq. (4) in the linear-response regime, we find that in the regions I and II the 
neural network behaves as a damped oscillator driven by a force F e (t), 

d 2 Ap e (t) „„ dAp e (t) ?1 . , 

£2 + 2 ^ C °°—^ + °^ A P^) = Pe{t), (10) 

(see, for example, in [18]). Here we introduced the damping ratio £ = y r /ft>o and a 
frequency C0q = (y 2 + yf ) l l 2 . In region I, the network is critically damped because £ = 1 
and it is underdamped in region II, where £ < 1. In the case f e (t) ^ and fi(t) = 0, 
the force F e (t) equals F e {t) = (1 — p e (°°))(B2f e (t) +df e (t)/dt). The parameter B2 was 
defined above. Solving Eq. (10) leads to a response function, 



X ee {t - = X.e-^-^ sin Wit - 1') + d> 



(ID 



where X e =(l-p e (oo))[l+(5 2 -7r) 2 /?f] 1 ^ 2 and 4> e =tan -1 [Yi/(B2—y r )] (one finds a sim- 
ilar result for X ; and 4> ; of inhibitory neurons). If y r > 0, then Eq. (11) shows loss of 
memory in the neural network with increasing time interval t—t'. If y r tends to zero, 
the memory becomes long-range. The Fourier transform % ee {oS) of the linear response 
function is 

= <^>^> . (12) 

(Oft- C0 2 + 2iQ(OqCO 

Equation (12) shows that at C, < 1, the neural network acts as a band-pass filter. The 
spectral intensity as a function of CO has a maximum at a resonance frequency <y r 
(O0a/1 — 2£ 2 at £ < 1/ \/2. The maximum value ||^ ee ((O r ) || 2 depends on the noise level 
(n). When approaching the critical point, y r — > 0, the value ||^ ee (w r )|| 2 diverges as 
% ee {(D r ) oc \/yj: — >■ oo, while the angular frequency of damped oscillations y tends to the 
frequency of stable network oscillations. This behavior signals that, in this regime, in the 
presence of noise, a neuronal network can amplify periodic signals. This amplification 
may be a mechanism of stochastic resonance observed in brain [19]. 

The band pass filter behavior described by Eq. (12) seems to be supported by mea- 
surements of response of rat CA1 neurons to afferent stimulation in vitro [12]. These 
measurements revealed that the majority of rat CA1 neurons act collectively like a band- 
pass filter and fire synchronously in response to a limited range of presynaptic firing 
rates (20 — 40 Hz) that are in the range of gamma oscillations in the rat hippocampus 
[20]. One can also note that, a long time ago, a number of characteristics of a band-pass 
filter behavior and a resonance response on sin wave trains already have been observed 
in EEG recordings of alpha activity [21]. Based on Eq. (12), we suggest that band-pass 
filter behavior observed in Sasaki et al. [12] and Tweel [21] is a manifestation of the 
critical phenomena near to the transition to neural network oscillations. 



STOCHASTIC FLUCTUATIONS OF NEURONAL ACTIVITY 

EEG measurements demonstrate that brain activity always contains a stochastic com- 
ponent. In this section we will show that stochastic fluctuations are enhanced when a 
neural network is close to the critical point of the non-equilibrium phase transition. For 
characterizing stochastic fluctuations, we introduce the autocorrelation function 

1 f T 

Cab(t) = - 8pa{h)8p b {h+t)dt U (13) 

1 Jo 

where 8p a (t) = p a (t) —p a describes fluctuations of activity p a (t) of population a, 
a = ej, around the mean value p a (see, for example, Ref. [22]). C a t,(t) is a measure 
of correlations between values of 8p a {t\) and 8pt(ti + 1) at two different instants 
separated by a lag t and averaged over an arbitrary large time window T. The Wiener- 
Khintchine theorem states that the power density spectrum of the fluctuations is the 
Fourier transform of the autocorrelation function. 

For calculating the autocorrelation function, one uses the standard method [22, 23]. 
In the deterministic equation (4), we assume that f a (t) is a stochastic force that satisfies 
conditions (f a (t)) = and (fa(t)fb(t')) = fo8(t — t')8 a ^. If fluctuations are small, the 



autocorrelation function may be found in the linear response theory [22, 23]. Assuming, 
for simplicity, fj(t) = 0, we obtain Eq. (9) that leads to 



C ee {t) = 2%fl 



\ e iC0t \\Xee(C0)\\ 2 dC0, 



(14) 



J — oo 



where the linear response function # ee (co) is given by Eq. (12). In the region of damped 
oscillations, the autocorrelation function C ee (t) has a form 



The parameter A e and the phase behave as A e oc \/y r and <J> e y r /y at small y r . 
For inhibitory neurons we obtain a similar behavior. Thus, stochastic fluctuations of 
activities of excitatory and inhibitory neural populations are enhanced when approaching 
the critical point y r = 0, Eq. (8), of the emergence of network oscillations (see Fig. 1). 
However, the linear-response approximation is not valid when fluctuations become 
sufficiently large. This occurs near to the non-equilibrium phase transition and non- 
perturbative methods are required for calculating C a b{t). 



In the present paper, using a cortical model with stochastic neurons, we have showed 
that, in neuronal networks, spontaneous appearance of sustained network oscillations 
occurs as a non-equilibrium phase transition. The critical point is determined by the level 
of noise, structure of the neural network, the balance between excitatory and inhibitory 
neurons, and other parameters. We have found critical and resonance phenomena that 
precede the transition. The important property of this transition is that, at the critical 
point, the relaxation time of the neuronal activity to a steady state becomes infinite in 
the infinite size limit. An increase of the response of neural networks to periodic afferent 
stimulations and a strong enhancement of stochastic fluctuations of activities of neural 
populations are also the critical phenomena that precede the transition. Note, that these 
phenomena are general properties of second-order phase transitions observed in physi- 
cal, chemical and biological systems (see, for example, Stanley [24], Haken [25], Kelso 
[8]). These critical phenomena have been observed near the non-equilibrium phase tran- 
sition in human hand movements [4, 5, 6, 7, 8]. The noise-induced nonequilibrium phase 
transition found in [26] is one more example of a phase transition with similar critical 
phenomena. Furthermore, we have demonstrated that near to the critical point, neuronal 
networks behave as damped harmonic oscillators or band-pass filters in agreement with 
band-pass filter behavior observed in vitro in networks of CA1 neurons in mammalian 
brain [12]. We suggest that band-pass filter behavior is a manifestation of critical phe- 
nomena near to the transition to network oscillations. 

We have also demonstrated that, in the cortical model, stochastic neural activity 
generated by a stochastic force is similar to spontaneous alpha activity observed in EEG 
recordings of both a normal man and human epileptic seizures of petit mal activity [27]. 




(15) 



CONCLUSION 
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